ESTIMATION OF ANISOTROPIC GAUSSIAN FIELDS THROUGH 

RADON TRANSFORM 
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Abstract. We estimate the anisotropic index of an anisotropic fractional Brow- 
nian field. For all directions, we give a convergent estimator of the value of the 
anisotropic index in this direction, based on generaHzed quadratic variations. We 
also prove a central limit theorem. First we present a result of identification that 
relies on the asymptotic behavior of the spectral density of a process. Then, we de- 
fine Radon transforms of the anisotropic fractional Brownian field and prove that 
these processes admit a spectral density satisfying the previous assumptions. Fi- 
nally we use simulated fields to test the proposed estimator in different anisotropic 
and isotropic cases. Results show that the estimator behaves similarly in all cases 
and is able to detect anisotropy quite accurately. 



Introduction 

The one dimensional fractional Brownian motion (IBm) was defined through a 
stochastic integral by Mandelbrot and Van Ness [20] in 1968 for the modeling of 
irregular data such as the level of water flows or economic series. Let us recall that 
this process is a Gaussian zero mean process with stationary increments character- 
ized by its so-called Hurst index H G (0, 1) and denoted by Bh = {BH{t);t G M}. 
A generalization of Bochner's Theorem allows to give a spectral representation of 
its covariance function, namely 

(1) Gov Bh{s)) = f (e-*« - 1) (e-« - l) ler'^-^e 

The function |^|~^^~^ is called the spectral density of the fBm. Processes with that 
kind of spectral density are called "1 //-noises" in the terminology of signal theory. 
The Hurst parameter is the index of irregularity of the fBm. It corresponds to the 
order of self-similarity of the process and to the critical Holder exponent of its paths. 
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In this paper we consider ci-parameter real-valued Gaussian fields with zero mean 
and stationary increments, defined through a spectral representation 

(2) |£ (e-^*-« - 1) fiO'^'Wmte , 

where ■ is the usual scalar product on M*^ and W is a complex Brownian measure 
with W{-A) = W{A) for any Borel set A C R'^. The function /, called the spectral 
density, is a positive even function of (M'^,min(l, d^). 

A natural extension of the 1-dimensional fBm is obtained when the spectral den- 
sity is given by l^l"^^""^, where | ■ | is the euclidean norm on M*^. This yields a zero 
mean Gaussian field with stationary increments which is isotropic and has therefore 
the same critical Holder exponent H in all directions of W^. 

In order to get an anisotropic field with stationary increments A. Bonami and 
A. Estrade define in [7] an anisotropic fractional Brownian field by considering a 
spectral density of the shape 

(3) A(o = \^r''^^^-', 

where the power h{^) E (0,1) is an even function which depends on the direction 
^/|.^| of M''. Other generalizations have been proposed for anisotropic data model- 
ing like the fractional Brownian sheet |19) or the multifractional Brownian motion 
introduced simultaneously in [5] and [23], where the Hurst parameter H is replaced 
by a function depending on the point t G M"'. However such generalizations yield 
models with non stationary increments. 

Here we consider an anisotropic fractional Brownian field X defined by ([2]) with 
([31) for some anisotropic index h and we focus on the identification of h. As al- 
ready noticed in [7j this index cannot be recovered by analysing X line by line 
since its regularity along a line does not depend on the direction. Hence, to recover 
anisotropy, authors give another directional analysis method, which is based on field 
projections. Actually, the critical Holder exponent of the process {RgX(t);t G M} 
obtained by averaging the field along an hyperplane orthogonal to a fixed direction 
9, called Radon transform of X, is proved to be equal to h{6) + for this direction. 
An estimator of h{9) is then proposed in [2j, using quadratic variations to estimate 
the regularity of the process RgX . However, no speed of convergence nor asymp- 
totic normality can be found under their weak conditions that the spectral density 
behaves like at high frequencies. 

Actually, many estimators for the Hurst parameter of a ID fBm have been pro- 
posed, based for example on time domain methods or spectral methods (see |9j and 
[3] and references therein). Quadratic variations can lead to relevant estimators 
with asymptotic normality of the Holder exponent of more general Gaussian pro- 
cesses with stationary increments as proved in [H] or [I7] for instance. Moreover 
in [I8] the authors give precise bounds of the bias of the variance and show that 
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minimax rates are achieved for this kind of estimators. However, these previous 
works need assumptions on the variogramme of the process X of the following type 

(4) v{t) = E ((X(t + to) - X{to)f) = + r{t) and r{t) = ^o^ , 

with further regularity assumptions on the reminder. This leads to a first restriction 
on the set of values of H since H must be in (0, 1). For this reason we adopt here 
a spectral point of view related to the problem of the identification of filtered white 
noise introduced in [4j. In the simplest case of this paper, the authors consider a 
spectral density / given by 

/(0 = c|e|-'^-' + i?(0, 

with G M+ \ N, c> and i? G C\R), satisfying \R^^\0\ < C\^\~^^^^~'~' , with 
s > 0, for < j < 2. It follows that, when H > 1, one has to consider a process no 
more with stationary increments but with higher order stationary increments. 

In this paper, under the assumption that the spectral density of a Gaussian process 
with stationary increments satisfies for some c, s > 

(5) /(o = cier'^"^ + , o d^r^^-^"^), 

we give estimators of H > using quadratic variations of the process. On the 
one hand, for H G (0,1) this assumption leads to ([4]), using the fact that v{t) = 
4 JjgSin^ (y) fiOd-^- On the other hand, H can be an integer and we do not need 
that the reminder is in C^(]R) as in [4]. 

Then, with further assumptions on the derivative of / at high frequencies, we get pre- 
cise error estimates and asymptotic normality. Our main result is that the spectral 
density of the Radon transform of an anisotropic fractional Brownian field satis- 
fies ([5]) with H = h{9) + Therefore we can get consistent estimators of the 
anisotropic index h with asymptotic normality using quadratic variations of the 
Radon transform process. 

The paper is organized as follows. The first section is devoted to the estimation 
of the Holder exponent of Gaussian processes {d = 1) with spectral density, using 
generalized quadratic variations. We give estimators with asymptotic normality un- 
der assumptions that rely on the asymptotic behavior of the spectral density. In the 
second part we estimate the anisotropic index of an anisotropic fractional Brownian 
field, using Radon transforms of the field. These transformations lead to Gauss- 
ian processes with spectral densities for which we give an asymptotic expansion. 
Then we apply the results of the first part to this process. In the third part, we 
test the proposed estimator using anisotropic and isotropic simulated fields in two 
dimensions. 
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1. Identification of the exponent for a ID-process 

We prove in this section a first identification result in a general setting. It is based 
on the well-known fact that a consistent estimator of the critical Holder exponent of 
a Gaussian process with stationary increments can be recovered using generalized 
quadratic variations (see [O] or [17] for instance). Actually, many authors have 
considered these estimators under assumptions based on the variogramme of the 
process. This is not adapted here for our framework and we prove similar results 
under assumptions based on the asymptotic behavior of the spectral density. Let 
us recall that, up to a constant, the spectral density of a ID fractional Brownian 
motion of Hurst parameter H G (0,1) is given by the function |^|~^^~^. Remark 
that no process with stationary increments can admit such spectral density when- 
ever H > 1 since this function does not belong to L^(]R, min(l, |.^p)(i.^) anymore. 
However one can obtain such spectral densities by considering processes with higher 
order stationary increments (see [5] and [24] for instance). Moreover there is no 
restriction to define a process X with spectral density asymptotically equivalent to 
that kind of function for any H > 0. Actually, when H > 1, writing H = j + s with 
J e N and s G [0, 1), X will be j times differentiable in mean square and X^^^ will 
admit s as critical Holder exponent. 

Let us consider a zero mean Gaussian process X = {X(t); t G M}, with stationary 
increments and spectral density / G L^(M, min (1, |^p)(i^). Let us assume that / 
satisfies (Ej) namely 

for some H,c,s> 0. We observe a realization of X at points for k = 0, . . . , N: 

Our purpose is to estimate H. The key idea of the former works concerning the 
estimation of Holder exponent of Gaussian processes with stationary increments 
as [m [171 [iHl E] for instance, is to consider increments of the process to get a 
stationary process. For instance, since X has stationary increments the process 
[X (^) — X (-^) ; t G M} is stationary. More generally, one can consider the fil- 
tered process of X 



ZnA^) = o,kX ( ) > for 

I — n \ / 



t G 



fc=0 

This defines a stationary process when a = (cq, . . . , a;) is a discrete filter of length 
I + 1 and of order K > 1 {l,K E N with I > K), which means that 

akk" = for < r < if - 1 and ^ auk^ ^ 0. 

fc=0 fc=0 
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For such a filter a, the observations of X allow us to compute Z^^aip) for p = 
0, . . . , N — I. Let us remark that if we choose the filter of order 1 given by a = (1,-1) 
then ZN,a{p) = X (^^) —X (-^) is the increment of X at point p/N with step 
More generally, for > 1, the increments of order i^' of X at point p/N with 
step are given by Ziy ^ip) for the filter a of order K and length K + 1 with 

ttk = {j^^ ^ k\{K-ky. for < A; < K. Let us point out that, since 

X has stationary increments, Zjva is a stationary process for any filter of order 
K > 1. Following [17] we also consider the filtered process of X with a dilated filter. 
More precisely, for an integer m > 1, the dilation a" of a is defined ior < k < lu 
by 

u _ j o,k' if k = k'u 
^ \ ^ otherwise. 

lu I 

Since ^Aj^a^ = u'^'^k^ak, the filter a" has the same order than a. Then, for 

fc=0 k=0 

a filter a of length / and of order K > 1 and m > 1, due to the stationarity of 
the corresponding filtered process Z^^ g^^, we can estimate the empirical variance of 
^Af,a"(0) based on the observations oi X by considering 

^ N-lu 

(6) Vn,.^ = ^ _ ^ ^ i^N,Ap)? > 

which we call generalized quadratic variations of X. 

Let us point out that we consider here the same kind of set of locations for the 
sum as in [171 E] but one could also consider more general one as done in [T8] . 
Moreover, let us remark that one can also consider m- variations of the process that 
estimate E (^Ar,a"(0)"^). We will focus here on the quadratic variations (m = 2). It 
is motivated by a result of J. F. Coeurjolly [10] who proves that, in the fractional 
Brownian motion case, the asymptotic variance of the Hurst parameter estimator is 
the lowest for m = 2. 



To build estimators of H the main idea is to choose a filter a such that 

2\ r^AT-2H„.2H 



u 



Then, by considering estimators given by 



(7) T, 



precise estimates of the asymptotic behavior of Gov (Zjv,a"(p)5 ^Af,a^(p')) allow to 
get the almost sure convergence of (T^r a", T^v.a^) to (1,1) with asymptotic normal- 
ity. Then, an asymptotic estimator of H can be built by considering for instance 
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V. 



The end of this section is devoted to the rigorous proofs of these statements, under 
assumptions that rely on the asymptotic behavior of the spectral density. 

Our first result allows to get an asymptotic development of E(VAf,a") according 
to the asymptotic development of the spectral density at high frequencies. Let us 
remark that if we associate to the filter a the polynomial 

I 

PJyX) = y^OfcX^, for X eR, 

then a is a filter of order K if and only if Pa\l) = 0, for < r < A' — 1 and 

Proposition 1.1. Let X = {X{t); t G M} 6e a zero mean Gaussian process, with sta- 
tionary increments and spectral density f . We assume that that there exists H > 0, 
c > and s > such that f satisfies ([5]). Then for any filter a of order K > 1 and 
u> 1 we have the following asymptotics: 

^2KfpQll\'j^2Kf^^^^^ + O (1) tfK<H 



^^^('^^) 2clogiV +^0^ (logiV) ifK = H. 



Moreover, if K>H then E^{H) = u^^ |Pa(e-^^«)| 1^1"^^"^^ < +oo and 

( O (Ar-2(i^-H)N ifK-H<s/2 
iV2^E(\/^,,u)=cn2^Ei(iJ)+i ^0_^(iV-MogiV) ifK-E = sl2 

"0°° (N-') ifK-H>s/2 

Proof. Since by assumption X has stationary increments and spectral density / it 
follows from ([2]) for (i = 1 that, when K 

ZNMt) [ e-^'^PJe-^^)f{0'^'W{dO. 

Therefore, for all p G Z, 

E = / Pa (e-^^^) ' f{i)di = E (V)v,a") . 
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Let US assume that K < H, since a is of order K, from Taylor formula, for e > 

|2 



there exists 6 > such that 



|Pa(e"*«)| 



l?l<^ 



P. ( e-^^ 



2K„,2K 



K\ 



< e^"^^. Hence, 



j\(\<^ 



According to ([5|), 

J\e\<^ 



Le^'fiOdC +^S-,J^^ if K<H 



and J, 



l?l>^ 



P. e~'it 



2c\ogN 



- O (1) iiK = H 



f{Cjdi= O (iV~2^), which yields the result. 

A''— »+oo 



Now, let us assume that K > H and write /(^) = c|^| ^ + R{C) with 

Re (R,min(l, satisfying = O . A change of 

|5|->+oo 

variables leads to 



P. e-'N 



with El{H) = J^\P^{e-'^)\ \^\~^"-^d^ < oo since H < K. Then, applying the 
previous results to R with H replaced by + s/2 we obtain the result for the 
reminder term, which concludes the proof. □ 

Therefore, to recover the parameter H , one has to consider a filter a of or- 
der K > H. In this case, for any u > 1, by stationarity of Z]^sc^ we obtain 
that E(V)v, N~'^^u^^cEl{H). In order to prove the almost sure con- 

vergence of (T/v,a"5 ^A^.a^) for > 1, with T/v,a" given by dTj), one has to estimate 
Gov (VAT^a", VN,a^). Actually, 

^ ^ N~luN~lv 

Gov (\^7V,a", ■^^W.av) = ^ ^ _ X] Gov (^7V,a"(p)^ ^7V,av(p')^) , 



AT 



p=0 p=0 



with Gov (Z^^,a"(p)^ ZAr,av(p')^) = 2G0V {Z n ,^,^{p) , Z N ,a:^{p')f siuCC (^Ar,a" (p) , ^A^.av (p' 

is a Gaussian vector. Moreover, 



Gov (^Ar,a"(p), ^JV,av(p')) 

Let us denote 



.■(p-p')s 



^Pa(e-^^)Pa(e-^^)/(e)(ie 



(8) r"^;(p) 



" ^ N 



f\i)di where = Pa (e'^^^) Pa (e^-^) 
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such that 

N~lu 
p=-N+lv 

It is obvious that under ([5]), for K > H we have 

N-lu 

However, we have to consider Yl ^Ifaip)"^ further assumptions have to be 

p=-N+lv 

done to see when {J^e-'P^hl^'' {^) \C\-^"-^d^) ^ is in f{Z). 



Proposition 1.2. Under the assumptions of Proposition M . 1\ if we assume moreover 
that f is dijjerentiable on M \ (— r, r), for r large enough and 

(9) /'(O = -i^H + 1)^^ + 0|,|^+^ {]^) ' 

then, for K > H and any 6 < min(2(A' — H),l) with 6 > max(l — 2if, 0), one can 
find C > such that, for all \p\ < N, 

KM<cN-'''{i + \p\r'. 

Moreover, for any K > H+l/A we have C^'" (H) = 2^ (4 e'^P^/i!;'^ (0 \^\~^^~^d^Y < +oo 
and 

Cov{Vn,.^,Vn,..) ~ c2cr(i7)iV-4^-i. 

A''— >+oo 

Proof Let K> H. Let us write /(O = c\^\~^^~^+R{0 with R e (M, min(l, \^\'^^)d^) . 
By a change of variables, we can write 

Let us focus on the first term and remark that for u = f , since /ia"(0 = \Pa (e~™^) |^ 
is real 

/ e-'^^hr{i)\i\-'"-'di= [ cos{pOh:^^{0\^\-'''-'d^. 
Jr Jr 

Moreover, since a is a filter of order K, by Taylor formula one can find C > such 

that 

d 



<Cmin(l,|er-i) 



(10) <Cmin(l,^2i^^ and 

Therefore, when p 7^ 0, we can integrate by parts 
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Then, for any 6 < min(2(A' - H),l) with 6 > max(l - 2H, 0), 



< \P\' 



1~S 



with 4 \^\^~^ ^ (/i^'" (0 \^\~^^~^) < +00 according to Writing h^' (0 as 
\[\P. (e-™«) + Pa (e"'^«) I' + ^ |Pa (e~*"«) + ^Pa (e"™«) I' - (1 + ^) (/^r + 



we also get that for any 5 < min(2(_R' — H), 1) with 5 > max(l — 2H, 0), for all 

peN, 

' <C(l + b|)-^. 



Jr 



Therefore {J^e ^P^h^'^ {^) \^\ ^'^Opez ^^(Z) as soon as one can choose 5 > 
1/2, which is possible when K > H + ^. It remains to study the second term. 
Let e > and choose r > such that R is differentiable on R \ (— r, r) with 
\R^^\0\ < e|Cr^^"^"^' for J e {0, 1} and |^| > r. We write 



f e-^'^hr(^)R{Od^= I + 



ICI>r- 



For the first integral, since R G L^(M, min(l, \E,\'^^)dS,), we remark that 



|C|<r 



< C{r)K 



-2K 



where C{r) is a constant that only depends on r, which may change line by line. 
For the second integral, for all p G Z, similar computations as previously leads to 



l?l>'- 



e ^^/^r ( ^ ) RiOd^ 



< C{r)N~^^ + Ce{l + \p\y'N' 



~2H 



for any 5 < min(2(/i — H), 1) with 6 > max(l — 2H, 0). Therefore, for any \p\ < N, 



Rim 

It follows that one can find C > such that, for all \p\ < N, 



2H 



~2H 



withN'^T-^lip) c 



Finally, when K > H + l/A one can choose 5 G (1/2,1) such that = 

2Epez [k^^'^^^hY (0 l^r^^^^^^O is fiiiit^- the dominated convergence theo- 
rem, 

N-lu 

iV^^+lC0v(\/;V,a",V^7V,av^ 



1 ' ' iv — ^+oo 



□ 
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We can now state our first identification result. 

Proposition 1.3. Let X = {X(t); t eM.} be a zero mean Gaussian process, with sta- 
tionary increments and spectral density f, which satisfies the assumptions of Propo- 
sition \1.1\ and Proposition Let a. be a filter of order K > H and u, f > 1 two 
integers with u ^ v. Then, almost surely, 

Hn,^.{u,v) = — — ^^log ( ) „^ — ^ H. 

2\og{u/v) \VN,e,^ J N^+co 

Moreover, for K > H + let us denote 



A\og\u/v)\El{HY El{Hy E-{H)El{H) 
Then, when s > ^, Vn (h^^{u, v) - h] ^ U (0, il'\H)) , with 

and, when s <\, 

Proof. Following classical computations on Gaussian quadratic forms as in |[14j for 
instance, from Proposition ll.il and 11.21 we obtain that for all i^' > if 

with for a\\K> H + 1/4, ViV {Tn,s^,Tn,^^) ^ M (0, T.''/'{H)) and 

' ' Af— »+oo 



iVCov((T;V,a",T;V,av)*(T^,a",r^,av)) ^ ^T\H) , 



where 



yu,v(rr. _ f Cl'^{H)/El{Hf Cr{H)/El{H)El{H) 
^a l^J ^ Cl^\H)/El{H)El{H) Cr{H)/El{Hf 

Then, using Taylor Formula for the function g{x,y) = log (see Theorem 3.3.11 
in Hi] for instance) we get that for K > H almost surely log ( ) — > with, 
for all i^: > if + 1/4, ViVlog (^^) AT (0, Tl'"{H)) and 





a" \ 




avy TV- 




^iV,a" 







2^ 



iVE(iog(;^) ] ^-^^^K^^'m 
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where T'^''"{H) = 41og(n/t;)^7^'^'(iJ). Then let us write fovuj^v 

HnAu, v) = — — -— - log — + log 

21og(n/i;) V VlE(VAr,av)/ V^A^.a- 

According to Proposition 11.11 it is straightforward to see that for K > H 

H]\faiu,v) — > Has. 

Moreover, ioi K > H + 1/4 and s > 1/2, Proposition ll.il leads to 



So in this case f i;) - iy) ^ A/" (0, 7^'^(i7)) and 

V ' / A''— »+oo 

Let us point out that if s < 1/2, from Proposition 11.11 using the fact that K > 
H +l/4> H + s/2 we get 



E(V)v,a") {u 



Therefore, 



E((if;raK.)-if)')=^P^^(iv--). 



□ 



This estimator is used in the next section in order to estimate the anisotropic 
index of an anisotropic fractional Brownian field. 

2. identification of the anisotropic index of anisotropic fractional 

Brownian fields 

In this section we consider an anisotropic fractional Brownian field 

X = {X{t) ; t G M^} , 

as introduced in [7], which is a zero mean Gaussian random field, with stationary 
increments and spectral representation ([2]). The spectral density is given by (JSj), 
namely fh{0 = where h is an even homogeneous function of degree 

with values in (0,1), called anisotropic index. To determine anisotropy of such a 
field one could try to estimate its directional regularity by extracting lines of the 
field along various directions. However, for anisotropic fractional Brownian fields, 
this method fails. Actually, when ^ is a fixed direction of the sphere S"^~\ one can 



12 



HERMINE BIERME AND FREDERIC RICHARD 



prove that the process {X(t9)]t G M} is still a zero mean Gaussian process with 
spectral density given by 

peR^ [ MpO + l)d7= f p'''^'^^'^^' (1 + |^|2)-M^+7)-<^/2 

where {6)-^ stands for the hyperplane orthogonal to 6. Therefore, according to 
Proposition 3.6 and Proposition 3.3 of [7] this process admits a critical Holder ex- 
ponent equals to the essential infimum of the function h. Then, the study of the 
generalized quadratic variations of such a process can at most allow us to recover 
ho. To deal with this obstruction and in order to study processes rather than fields, 
the authors of |7j have introduced the Radon transform of such fields. 
When a function / is integrable over M"^, one can define its Radon transform on M 
(see [27] for instance), in the direction 6^, by 

7^e/(^) = f f{s + t9)ds, for all t E R. 

For a function / which does not decay sufficiently at infinity, one can integrate it 
against a window. Let p be a smooth function defined on {6)-^, that compensates 
the behavior at infinity of /. Then one can define the windowed Radon transform 
of / on M, in the direction 6^, by 

nepf{t)= I f{s + te)p{s)ds, iov a\\t eR. 

We should use strong assumptions on the anisotropic fractional Brownian field X in 
order to define its windowed Radon transform as an integral. However, according to 
Proposition 4.1 of [7], one can define the Radon transform of X, with a convenient 
window p, in the direction 6, by a discretization of the integral. For notational sake 
of simplicity we deal with the direction 6* to be = (0, ... , 0, 1) and identify the 
space W^^^ X {0} to R'^"^. Let us choose p a function of the Schwartz class S (M'^"^), 
with real values, ie p is a smooth function rapidly decreasing 

(11) ViV e N, Vx e M"'-\ |p(x)| < C^(l + |x|)-^. 

Then, the process 

2-„(d-i) X(s,t)p(s), Vt G M, 

admits a limit in L^(f2) for the finite dimensional distributions, when n tends to 
infinity. This limit is called the Radon transform of X with the window p and is 
denoted by RpX = {RpX{t); teR}. 

Let us remark that one can define the Radon transform of X, with the window 
p, under less restrictive assumptions on p, as soon as the previous limit exists. The 
existence of the limit process is proved in both and [6] and relies on the slow 
increase of the covariance function of X due to its stationary increments and on its 
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mean square continuity. 

Let us also point out that one can define the Radon transform of X with the window 
p for any direction 6. Actually, it is sufficient to chose kq a rotation of that maps 
6*0 = (0, . . . , 0, 1) onto 6. Since X o kq i?, still an anisotropic fractional Brownian 
field, with anisotropic index given by ho which satisfies the same assumption as 
h, we just have to consider the Radon transform of this field with the window p. 

By linearity of such a transformation, the Radon transform of X with the window 
p is still a zero mean Gaussian process with stationary increments and it admits a 
spectral density given by the Radon transform of the spectral density fh of X, given 
by (l3|, against the window \p\^ , 

(12) 7^|^./,(p)= / /,(7,p)|p(7)|2d7, for allpGM, 

where p is the {d — l)-dimensional Fourier transform of the window p. To estimate 
the Holder regularity of this process, we will use its generalized quadratic variations 
as introduced in the section [H Therefore, we have to study the asymptotic behavior 
of 7^l^2//i in order to apply Proposition II .Si We prove and use the following general 
result on the windowed Radon transform. 

Proposition 2.1. Let h and c be given functions on M"'. Let a > 0. We assume 
that h and c are even homogeneous functions of degree 0, Lipschitz of order a on 
the sphere, with h positive. 

Let 6o > 0. Let f be a function defined on M.'^ such that, for all 6 G (0, 6q), 

Choose p E S (M^^^^ such that J^d-i p{l)d'y = 1. Then, the Radon transform of f 
with the window p satisfies, for all 6 G (0, 6i), 

T^pfip) = + o (^p^^) when peRand \p\ +oo, 

with 6i = min (5o, a) ■ 

Proof. Let p be a function of S (M"'"-'^) with J^d-i p{l)d'y = 1. For p G M, with \p\ 
large enough, one can define the integral 

T^pfip) = / f{l^p)p{l)dl- 

We want to estimate its asymptotics when \p\ +oo. First, let us assume that 
there exists A> 1 such that, for ^ eW^ and |^| > A, 

/«) = ^ 
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with h and c satisfying assumptions of Proposition I2.1[ In this case, we will prove 
that for all < (5 < a, 

(13) Tlpfip) = fipOo) + odpr'^^^")-^) when p G M and \p\ +oo. 

For \p\ > A, since /jgd-i p{l)dj = 1, let us write 

n,f{p) = f{p9o) + [ Uil.P) - fipdo)) p{l)di. 



Then, it is enough to give an upper bound for 

{f{i.p)-f{pe^))p{i)di. 



Let us denote h {^S'^ ^) = [/iq, hi] with /iq > by assumption on h. Since p is rapidly 
decreasing, for all s > and G N, 

(/(7,P) - f{pdo)) p{l)dl = 0|,|^+oo(br'«"^1, 

7l>|ph 



which is negligible compared to \p\ ^'^^"^ ^ as soon as iV > 
Thus, it is sufficient to consider 



s 



A.(p)= / (/(7,p)-/(p^o))p(7)rf7- 

l7l<lph 



But, 



|A,(p)| < / |c(7,p)| 

|7|<H'' 

1 



Q^|2 _^p2)^(7,P)/2) |p|Meo) 

c(7,P) -c(^o)||p(7)M7- 



|p(7)M7 



IP^'^ Jh\<\P 
Let us use the Lipschitz assumptions on h and c. 



Lemma 2.2. //(? is an homogeneous function of degree 0, Lipschitz of order a on 
the sphere S'^^^, then there exists C > such that for all p ^ Q and 7 G M'^"^, 

lfi'(7,P) - S-lO,!?)! < Cmin ( ( 14 ) 



bl 

Proof. The function g is continuous on the sphere and thus it is bounded. Then, for 

p ^ and 7 G M'^~\ 

|^(7,P) -^(0,p)| < 2||^||oo. 
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Moreover, g is Lipschitz of order a on the sphere. Then, there exists C > such 
that, for p 7^ 0, 

(7,p) (0,p) 



I^(7,P) -^(0,p)| < 



But, 



(0,p) 



2\ -1/2 



(|^|2+p2)l/2 

which concludes the proof of Lemma [221 



2 1 



< 



p2 



□ 



Let us recall that c is an even homogeneous function thus c(6'o) = c(0,p), for all 
p 7^ 0. Then, since c is Lipschitz of order a, one can find Ci > such that 



As,2(p) 



\V 



hieo) 



:(7,p) - c(0,p)||p(7)M7 < Cibr^(^«)-"(^-^), 



'y|<|p|'i 



which is negligible compared to \p\ ^ as soon as 5 < «(! — s) 
It remains to consider 



|c(7,p)l 



7|<|p|' 
1 



1 



bl 



hido) 



(|7|2 + p2)M7,p)/2 |^|Meo) 
I IM^^O) 



|7|<|P|' 



Let us write 
where, for p 7^ 0, 



\p\ 



h{Bo) 



(|^|2+p2)M7,p)/2 



|7r +p 



Ap) 



|p(7)M7 
\p{l)\dl- 



M^o)lnb|--M7,p)ln(bP+|7r) 
{h{0, p) - /i(7, p)) In IpI - ^/i(7, p) In ( 1 



I7r 



writing /i(6'o) = h{0,p) since /i is an even homogeneous function. Since h is Lipschitz 
of order a, by Lemma 12^21 for s < 1, there exists C2 > such that, for \p\ > A > e 
and I7I < IpI"*, 

In \p\ + 



< 2C2 



m 



Inbl < 2C2\p\ 



-a{l-s) 



In 
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The function t i— > In |t| tends to at infinity. Thus one can find Ag > 

such that for \p\ > As we get \l{p)\ < 1. Then, for \p\ > Ag 

|e^(p) _ i| < e|/(j9)| < 2eC2br"^'"'^ In bl, 

and finally 

A,,i(p) < 2eC2||c|Ubr'^'°^-"^'-^)ln \p\. 
For 6 < a, since \ As{p)\ < ^^sAp) + A^,2(p), taking s G (0, ^) C (0, 1) we get 

1 



As{p) = 0|p|. 



+ + 00 



and (fT3l) follows. In the general case, let us assume that, for all 6 G (0, 6q) and 

Replacing p by |p| and hhy h + 6 in the special case above, we get the result for the 
remainder. □ 

Let us remark that to simplify the statement of this proposition, we assumed the 
window to be in the Schwartz class. However it is proved in [6] p. 85 that the result 
still holds for a window p G L^(]R^"^) that satisfies 

(14) |p(7)| = O (p^i^) when 7 G M^"^ and I7I ^ +00, 

with M > hi-ho for h {S'^^^) = [ho, h]. In this case Si = min {5o, a^^^^f^). 



We can now state our main result concerning the identification of the anisotropic 
index of an anisotropic fractional Brownian field. We keep the notations of part [T] 
for the generalized quadratic variations of a ID-process and recall that we fix the 
direction ^0 = (0, . . . , 0, 1) G S'^~\ 

Theorem 2.3. Let X = {X(t);t G M'^} be an anisotropic fractional Brownian field, 
with anisotropic index h given by an even homogeneous function of degree with 
values in (0, 1), which is assumed to be in [S'^'^) . 

Let p be a window in S (IR*^"^) . Let RpX be the Radon transform of the field X with 
the window p. 

Let a be a filter of order K and u,v >1 two integers with u ^ v. If K > h{9o) + 
then almost surely 

= :^i4TT'°g ( v"fRx\ ) - ^ 

21og(u/t;) \VN,a-{RpX) J 2 TV-^+oo 
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Moreover, if K > h{9o) + ^ + 1/4 



with 



N (h^o) {u, v) - h{e,)^ AT ( 0, 7r ( He,) + 



d-i 

2 



Proof. It is sufficient to prove that the spectral density of RpX satisfies the assump- 
tion of Proposition II .Si From (fT2l) . this spectral density is given by the function 

n,^^fhip)= ! fh{{l.p))\p{l)?d^. for all PGM, 



with fh given by (l3|). Since we can divide RpX by a constant, we can assume that 
|p(7)|^'^7 = 1- Then, since 2h + d satisfies assumptions of Proposition [2TT] with 
a = 1, by Proposition 12.11 for all 5 < 1, 



We can also write this as 

T^\p\^fh{p) = |^|2(h,(eo)+^)+i + IpH+oo \ \p\2(hieo)+^)+i+s ^ 

Moreover, since h E (5^^"^), the function (7(7) := h{'y, 1) is differentiable on M'^"^ 
and, for 7,x G M''"^, 

1 / 7 • X 

(171"= + 1)^/"= KtHTT \ |7|^ + 1 

such that 

\D^g{x)\ < 2\\Dh\\^\x\. 

Thus, the spectral density fh is differentiable on W^' \ {0}. Let (7,^) G W^^^ x 
since h{^,p) = g{'y/p), we get 

p{2h{-f,p) + dy 



|/.(7,P) = Ml,p) (-^) In (|7P +/) - 



|7p + p"^ 



It follows that the spectral density 'R-^^-ifh is differentiable on M \ {0} and for p 7^ 
we have 

(7J„./,.)'(rt=£_,|A(7,p)l?(7)r^7. 

Let us write 

d I7I 

'O^fhil^p) = -^Fi(7,p) -pF2(7,p), 
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with 



and 

Therefore, 



F2(7,p) = (2/i((7,p)) + rf)A+i(7,p). 



p^ 

Then, by Proposition 12.11 whenever 6 < 1, 

2h{eo) +d (I 



T^\o\^F2{p) 



Moreover, for any e > small enough 

|Fi(7,p)| <2||D/i|UA-.(7,p)- 

Since I7I |p|^ is integrable over M'^"^ and rapidly decreasing, following the same lines 
as in the proof of Proposition 12.11 we get 

This allows us to conclude that 



^1^2 a) (p) 



2h{9o) + d 



d-1 
2 



Therefore RpX satisfies the assumption of Proposition 11.31 with H = h{9o) + 
and s < 1, which concludes the proof. □ 

Let us remark that if we choose a window such that \pf only satisfies (fTlIl for M > 

2{hi — ho) with /i (S''^~^) = [/io,/ii], the estimator hiy^a{(^o){u,v) still tends almost 
surely to h{9o). However the speed of convergence will depend on M. Actually the 
result of Theorem O holds when Af > l + 4(/ii-/io) whereas for M < 1 +4(/ii -/iq), 
Proposition 11.31 shows that, for all s < ^^^—jfp^^, 



N^(hZ(So){u,v)-h{eo)) 



Let us point out that we have restricted h to have values in (0, 1) so that the 
anisotropic fractional Brownian field is well defined. This can be weakened by 
considering fields with higher order stationary increments or with spectral density 
asymptotically of the order of fh- Similar results can be obtained with the further 
assumption that the partial derivatives of order 1 are asymptotically of the order of 
the partial derivatives of fh (see [6] for example) . 
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For numerical applications, one has to approximate the Radon transform of X. 
Following [2] we replace RpX{t) by 

(15) JmW = M-^'^-^) X{M-h,t)p{M-h), 

for p a smooth window with compact support and M > 1 an integer. Let us denote 
for a a filter of order K and length / and < M , 



(16) 



p=0 \fc=0 



The key point of the proof is to estimate the error due to the approximation of 
VN,a{Rp^) by Tlf^{X). Let us denote H = h{9o) + Under the assumptions 
of Theorem 12.31 following the same lines as in [2], for a > with a < ho = min/i, 

one can prove that there exists a positive finite random variable C such that, for all 
iV > 1, 



1/2 



< CM-" a.s. 



Moreover, since for K > H by Proposition 11.11 ; -y — ' 

find a positive finite random variable C such that 

V/v,a {RpXy^^^ < C'N" a.s. 
Then, for M'^N^ small enough, writing 

r^^jx)V2_v]v,a(/?pX)^/^ 



1 a.s., one can 



log 



2 log 1 



1/2 



one can find a positive finite random variable C" such that, a.s. 



log 



We can state the following result. 

Proposition 2.4. We keep the assumptions of Theorem \2.3\ and take p with compact 

min h 



support. Let H = h{9o) + ^^ and ho = min /i. If K > H and q > H/ho then, almost 



surely, 



2 log(n/i;^ 



log 



T^:a4X) \ d-1 



2 N^+c 



h{do) 



Moreover, if K > H + 1/4 and q > {H + 1/2) /ho then 



NlHN^M{u.v)-h{eo) 
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3. Numerical Study 

In this section we present a preliminary evaluation of projection-based estimators 
studied above. We first describe the synthetic datasets used for our evaluation and 
then discuss our estimation results. 

3.1. Simulation. A lot of numerical methods have been proposed these last years 
to simulate 1-dimensional fractional Brownian motion (fBm). Most of them give rise 
to approximate syntheses, such as the midpoint displacement method (see [22], for 
instance), the wavelet based decomposition ([21], [1], ^26j, etc.), or more recently 
a method based on correlated random walks [13]. A few of them can be applied 
not only for 1-dimensional fBm but also to simulate 2-dimensional (anisotropic) 
fractional Brownian fields and lead to approximate syntheses. Of course there exist 
exact synthesis methods based on the Choleski decomposition of the covariance 
function. This yields numerical problems due to the size of the matrix. In order to 
have fast synthesis one can use the stationarity of the increments by applying the 
embedding circulant matrix method [12j. By this way, we easily obtain fast and 
exact synthesis of 1-dimensional fBm t25j. Some authors, as in [8] and [16], apply 
this method for higher dimension but this does not yield to exact synthesis. Finally, 
M. L. Stein proposed a fast and exact synthesis method for isotropic fBm surfaces 
in [28]. 

For a preliminary evaluation, we used this method for generating a dataset con- 
taining three subsets of 1000 fBm fields simulated with each of the Hurst parameter 
values: h = 0.2 (low regularity) h = 0.5 (medium regularity) h = 0.7 (high regu- 
larity). fBm fields were generated on a discrete grid G of size (M + 1) x (M + 1) 
{M = 2^ = 512) which is defined by 

(17) G=|(|,|),o<*,.fa<M}. 

The method was implemented in collaboration with A. Fraysse and C. Lacaux, and 
the corresponding matlab codes are available at |ht tp : / / ciel . ccsd . cnrs . f r 

However, this method can not be generalized to simulate anisotropic fields. Hence 
we also used a spectral representation approximation (SRA) technique to generate 
anisotropic fractional Brownian fields (afB). The regularity of two-dimensional afB 
we generated differs in both vertical and horizontal directions. The spectral densities 
of these fields are of the form 

(18) v« = (6,6) € . {(0,0)}. no = { ii^^i<^i^^ 

In this expression, hh and /i„ form a pair of parameters in (0, 1) which characterize 
the anisotropy of generated fields. Their regularity is then given by mm{h^, hh). 
We used an approximation of the spectral representation in Equation ([2]) which was 
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obtained by discretization. The discrete approximation is given by 
with 

(19) y(^.|)=^ E E z{n,,n,)g(nn,,7rnX-^(-^'^^-^'^\ 

^ ^ ni=~M+ln2=-M+l ^ ^ 

where z{ni,n2) are (2M)^ independant realizations of complex random variables 
whose real and imaginary components are two uncorrelated zero-mean standard 
gaussian variables and the function g{x, y) = /^/^(x, y) for all {x, y) \ {(0, 0)} 
and g{0, 0) = 0. From a practical point of view, the sum in Equation (fTOl ) can 
be interpreted as a filtering in the Fourier domain of a white noise 2; by a low-pass 
linear filter characterized by the transfer function g. Based on this method, afB fields 
approximations can be easily and quickly simulated using the fast Fourier transform. 
Similar simulations are used in |15| for an evaluation of a different estimator. Such 
an approach can also be extended to the simulation of 3D fields and anisotropic 
fields having more than two different directional regularities. As an illustration, 
some simulation examples are shown on Figure [TJ 




(a) (b) (c) 

Figure 1. Some simulations of afB surfaces. Horizontal and vertical 
regularity parameters used for simulations are 0.7 and 0.7 in (a), 0.7 
and 0.5 in (b) and 0.7 and 0.2 in (c). 



By this way, we generated a second evaluation dataset containing six subsets of 
1000 fields on the grid G for each of the following regularity parameter pairs: low 
isotropic regularity {hh = 0.2, = 0.2), medium isotropic regularity {hh = 0.5, hy = 
0.5), high isotropic regularity {hh = 0.7, = 0.7), high/medium anisotropic reg- 
ularity {hfi = 0.7, = 0.5), high/low anisotropic regularity [hh = 0.7, hy = 0.2), 
medium/low anisotropic regularity {hh = 0.5, hy = 0.2). 
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3.2. Estimation. Given a field realization x, we first computed discrete Radon 
transforms Ih and in both horizontal and vertical directions, as defined in Equation 
(fTSj) for a window function of the shape l[o,i]: for all < /c, / < M, 




Then we computed generalized quadratic variations T^i for both directions (e = 
V, h), as defined in Equation (fT6l) . with the second order filter = (1, —2, 1) and a 
step of size 1/N where = M/2'^, for in {0, 1, 2, 3}: 

We also computed generalized quadratic variations T^2 with the dilated second order 
filter = (1,0,-2,0,1): 

M/2'^—4: 2 

Finally, we obtained the projection-based estimate h'^ of the index he of X in the 
direction e {e = v, h) for G {0, ■ ■ ■ , 3} as 

^-24)'- (1)4- 

For the evaluation, we computed empirical biases and standard deviations of es- 
timators over sets of simulated fields having same characteristics. Empirical biases 
were obtained as a difference between the mean parameter estimates and the real 
parameter value. In order to enhance the estimator ability to capture anisotropic 
properties, we also computed biases and standard deviations of differences between 
horizontal and vertical estimates. Whole results are reported in Tables [T] and [2l 

On Table [U we observe that standard deviations of estimation errors increase 
as the subsampling factor 2^ increase, meaning that subsambling of the projected 
signal has an effect on the estimator stability. However, standard deviations do not 
vary significantly when the subsampling factor is fixed and parameter values are 
changed. For instance, when z/ = 0, standard deviations are 0.049 and 0.055 when 
parameters values are 0.7 and 0.2, respectively. The order of standard deviation 
variations is about 10"^, for any fixed subsampling factor. 

Estimators underestimate the real parameter value. The underestimation in- 
creases as the parameter value decreases. For instance, estimation biases obtained 
for 1/ = are —0.047, —0.092 and —0.239 when parameter values are 0.7, 0.5 and 
0.2, respectively. Subsampling the projected signal reduces the underestimation bias. 
For instance, when h = 0.2, the bias is reduced to —0.041 when the subsampling 
factor is increased to 2^. 
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h 


7 / 






hi -\—/Ti 


0.7 





-0.047±0.049 


-0.045±0.045 


-0.002±0.069 


0.7 


1 


-0.012±0.056 


-0.017±0.061 


0.005±0.085 


0.7 


2 


-0.03±0.081 


-0.023±0.093 


-0.007±0.124 


U. I 


Q 
O 


n nc^/i -i-fi 1 1 Q 
-U.UJ34±U.iio 


n nA-i-Ct 11/1 
-U.U4±U.ii4 


n ni /i-un i kq 
-U.Ui4±U.ioo 


0.5 





-0.092±0.052 


-0.095±0.052 


0.003±0.073 


0.5 


1 


-0.034±0.069 


-0.035±0.071 


0.001±0.099 


0.5 


2 


-0.007±0.099 


-0.031±0.095 


0.024±0.132 


0.5 


3 


-0.029±0.131 


-0.023±0.136 


-0.006±0.192 


0.2 





-0.239±0.055 


-0.245±0.062 


0.006±0.082 


0.2 


1 


-0.112±0.079 


-0.131±0.082 


0.019±0.119 


0.2 


2 


-0.041±0.113 


-0.039±0.12 


-0.002±0.162 


0.2 


3 


0.002±0.163 


-0.039±0.147 


0.041±0.238 



Table 1 . Evaluation of directional regularity estimators on synthetic 
fBm surfaces of Hurst index h simulated using the Stein's method of 
exact synthesis. Values bh and b^, are empirical biases of horizontal and 
vertical regularity estimators and values ah and are their associated 
standard deviations. Values „ = bh — b^ are differences between bias. 
Values (j/i „ are standard deviations of differences between horizontal 
and vertical regularity estimates. 



On Table we observe that estimation biases obtained on fields simulated using 
the SRA technique differ from those of Table [U Indeed, these biases indicate that 
estimators overestimate parameter values. This is probably due to the SRA tech- 
nique which generates fields which are smoother than what they should be. The 
best evaluation of estimator biases are those of Table [T] which were obtained on 
exact simulations of fBm. However, we can notice that standard deviations of Table 
[2] are consistent with those of Table [TJ Despite simulation errors of the SRA tech- 
nique, we can rely on results of Table [2] to get a sense of the estimator properties 
when fields are anisotropic. In particular, we see that estimator biases and stan- 
dard deviations do not vary significantly from isotropic to anisotropic simulations. 
When the subsampling factor is fixed, standard deviations are about the same for 
isotropic and anisotropic simulations. On anisotropic simulations, biases still vary 
when estimated parameter values are changed. But, the bias of a regularity esti- 
mate in one direction does not depend on the parameter value estimated in the other 
direction. For instance, when z/ = 0, biases for the estimation of parameter value 
h = 0.2 are about 0.15 in all the simulation cases involving an index value of 0.2 
{{hhyhy) = (0.2, 0.2), (0.7, 0.2), (0.5, 0.2)). In conclusion, the field anisotropy seems 
not to have any effects on the estimator stability. 

Besides, we observe on the last column of Tables [U and [2] that standard deviations 
of estimate differences are about the same when u is fixed. In particular, when = 0, 
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standard deviations are about 0.065. On isotropic cases, biases are about 0. This 
suggests that it is possible to distinguish between isotropic fields and anisotropic 
fields for which absolute differences between horizontal and vertical regularities are 
above 0.065. 



7 

hh 


7 

K 


u 


7 1 

Oh±o-h 


7 1 

Ov±0'v 


7 1^ 

0h,v±O-h,v 


0.7 


0.7 





0.068±0.041 


0.069±0.041 


-0.001±0.06 


0.7 


0.7 


1 


0.003±0.063 


0. ±0.059 


0.003±0.087 


0.7 


0.7 


2 


-0.012±0.09 


-0.014±0.087 


0.002±0.126 


0.7 


0.7 


3 


-0.021±0.125 


-0.024±0.13 


0.003±0.182 


0.5 


0.5 





0.1±0.046 


0.102±0.044 


-0.002±0.065 


0.5 


0.5 


1 


0.012±0.07 


0.009±0.067 


0.003±0.097 


0.5 


0.5 


2 


-0.007±0.1 


-0.008±0.097 


0.001±0.139 


0.5 


0.5 


3 


-0.013±0.142 


-0.015±0.146 


0.002±0.207 


0.2 


0.2 





0.156±0.052 


0.16±0.05 


-0.004±0.073 


0.2 


0.2 


1 


0.034±0.08 


0.03±0.078 


0.004±0.112 


0.2 


0.2 


2 


0.004±0.113 


0.004±0.112 


0.±0.158 


0.2 


0.2 


3 


-0.004±0.163 


-0.007±0.164 


0.003±0.238 


0.7 


0.5 





0.071±0.041 


0.1±0.044 


-0.029±0.061 


0.7 


0.5 


1 


0.001±0.064 


0.002±0.067 


-0.001±0.095 


0.7 


0.5 


2 


-0.015±0.089 


-0.005±0.101 


-0.01±0.133 


0.7 


0.5 


3 


-0.026±0.131 


-0.014±0.137 


-0.012±0.189 


0.7 


0.2 





0.072±0.041 


0.157±0.052 


-0.085±0.065 


0.7 


0.2 


1 


-0.002±0.061 


0.029±0.078 


-0.031±0.1 


0.7 


0.2 


2 


-0.014±0.087 


0.01±0.114 


-0.024±0.14 


0.7 


0.2 


3 


-0.022±0.128 


-0.009±0.163 


-0.013±0.21 


0.5 


0.2 





0.098±0.045 


0.159±0.053 


-0.061±0.069 


0.5 


0.2 


1 


0.006±0.072 


0.032±0.079 


-0.026±0.108 


0.5 


0.2 


2 


-0.003±0.103 


0.007±0.117 


-0.01±0.16 


0.5 


0.2 


3 


-0.002±0.142 


-0.009±0.163 


0.007±0.211 



Table 2. Evaluation of directional regularity estimators on synthetic 
afB surfaces simulated using the SRA method of approximate synthe- 
sis. Values hh and h^ are horizontal and vertical regularity parameters 
used for field simulations. Values bh and 6„ are empirical biases of 
horizontal and vertical regularity estimators and values ct/j and 0"^, are 
their associated standard deviations. Values bh^v = bh — b^ are differ- 
ences between bias. Values ah,v are standard deviations of differences 
between horizontal and vertical regularity estimates. 
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